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Abstract 

Since its emergence in 1968, influenza A (H3N2) has evolved extensively in 
genotype and antigenic phenotype. Antigenic evolution occurs in the context 
of a two-dimensional 'antigenic map', while genetic evolution shows a char- 
acteristic ladder-like genealogical tree. Here, we use a large-scale individual- 
based model to show that evolution in a Euclidean antigenic space provides 
a remarkable correspondence between model behavior and the epidemiologi- 
cal, antigenic, genealogical and geographic patterns observed in influenza virus. 
We find that evolution away from existing human immunity results in rapid 
population turnover in the influenza virus and that this population turnover 
occurs primarily along a single antigenic axis. Thus, selective dynamics induce 
a canalized evolutionary trajectory, in which the evolutionary fate of the in- 
fluenza population is surprisingly repeatable and hence, in theory, predictable. 
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Author Summary 



Each year, hundreds of milhons become sick with the influenza virus and hundreds of 
thousands die from the infection. After recovery from the flu, a person gains permanent 
immunity speciflc to the infecting strain. However, owing to its RNA makeup, mutations 
occur rapidly to the virus genome. Some of these mutations change the shape of proteins 
visible to the human immune system, and thus alter the virus's antigenic phenotype. These 
mutations allow the virus to re-infect those who have previously recovered from an earlier 
strain, and thus will quickly spread through the virus population. It is because of in- 
fluenza's rapid antigenic evolution that the flu vaccine needs frequent updating. However, 
despite strong pressure to evolve away from human immunity and to diversify in antigenic 
phenotype, influenza, and especially influenza A (H3N2), shows paradoxically limited ge- 
netic and antigenic diversity present at any one time. Here, we propose a simple model 
of influenza that displays rapid evolution but low standing diversity and simultaneously 
accounts for the epidemiological, genetic, antigenic and geographical patterns displayed by 
the virus. We flnd that in evolving directly away from existing human immunity, the virus 
severely limits its future evolutionary potential. 



Introduction 

Epidemic influenza is responsible for between 250,000 and 500,000 global deaths annually, 
with influenza A (and in particular, A/H3N2) having caused the bulk of human mortality 
and morbidity [I]. Influenza A (H3N2) has continually circulated within the human popu- 
lation since its introduction in 1968, exhibiting recurrent seasonal epidemics in temperate 
regions and less periodic transmission in the tropics. During this time, H3N2 influenza has 
continually evolved both genetically and antigenically. Most antigenic drift is thought to be 
driven by changes to epitopes in the hemagglutinin (HA) protein [2] . Phylogenetic analysis 
of the genetic relationships among HA sequences has revealed a distinctive genealogical 
tree showing a single predominant trunk lineage and side branches that persist for only 1-5 
years before going extinct [3]. This tree shape is indicative of serial replacement of strains 
over time; H3N2 influenza shows rapid evolution, but low standing genetic diversity. 

This observation has remained puzzling from an epidemiological standpoint. Antigenic 
evolution occurs rapidly and strong diversifying selection exists to escape from human im- 
munity; why then do we see serial replacement of strains rather than continual genetic and 
antigenic diversification? Indeed, simple epidemiological models show explosive diversity 
of genotype and phenotype over time [H |5] . Previous work has sought model-based ex- 
planations of the limited diversity of influenza, relying on short-lived strain-transcending 
immunity [U [5] , complex genotype-to-phenotype maps [6] or a limited repertoire of anti- 
genic phenotypes [7]. 

Experimental characterization of antigenic phenotype is possible through the hemagglu- 
tination inhibition (HI) assay, which measures the cross-reactivity of HA from one virus 
strain to serum raised against another strain j8]. The results of many HI assays can be 
combined to yield a two-dimensional map, representing antigenic similarity and distance 
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between strains as an easily visualized and quantified measure \9\ . The path traced across 
this map by influenza A (H3N2) from 1968 until present is largely linear, showing serial re- 
placement of one strain by another; there are no major bifurcations of antigenic phenotype 

Herein, we seek to simultaneously model the genetic and antigenic evolution of the influenza 
virus. We represent antigenic phenotypes as points in a A^-dimensional Euclidean space. 
Based on the finding that a two-dimensional map adequately explains observed antigenic 
distance between strains [9j, we begin with antigenic phenotypes as points on a plane, but 
relax this assumption later on in the analysis. After exposure to a virus, a host's risk of 
infection is proportional to the Euclidean distance between the infecting phenotype and the 
closest phenotype in the host's immune history. Mutations perturb antigenic phenotype, 
moving phenotype in a random radial direction and for a randomly distributed distance. We 
implemented this geometrical model in a large-scale individual-based simulation intended to 
directly model the antigenic map and genealogical tree of the global influenza population. 
The simulation includes multiple host populations with different seasonal forcing, hosts 
with complete immune histories of infection, and viruses with antigenic phenotypes. As 
the simulation proceeds, infections are tracked and a complete genealogy connecting virus 
samples is constructed. Results shown here are for a single representative simulation of 40 
years of virus evolution in a population of 90 million hosts. 



Results 

The virus persists over the course of the 40- year simulation, infecting a significant fraction 
of the host population through annual winter epidemics in temperate regions and through 



less periodic epidemics in the tropics (Figure 1 A.). Across replicate simulations, we observe 



average yearly attack rates of 6.8% in temperate regions and rates of 7.1% in the tropics, 
comparable with estimated attack rates of infiuenza A (H3N2) of 3-8% per year \10\ 111] . 
Over the course of the simulation, the virus population evolves in antigenic phenotype 
exhibiting, at any point, a handful of highly abundant phenotypes sampled repeatedly 
and a large number of phenotypes appearing at low abundance ( Figure T^ ). By including 



measurement noise on antigenic locations, we approximate an experimental antigenic map 
of H3N2 influenza ( Figure Tp ). The appearance of clusters in the antigenic map comes 



from the regular spacing of high abundance phenotypes combined with measurement noise. 
Over time, clusters of antigenically similar strains are replaced by novel clusters of more 



advanced strains (Figure SI A.). Across replicate simulations, clusters persist for an average 
of 5.0 years measured as the time it takes for a new cluster to reach 10% frequency, peak 
and decline to 10% frequency. The transition between clusters occurs quickly, taking an 
average of 1.8 years. 

Remarkably, although antigenic phenotype is free to mutate in any direction in the two- 
dimensional space, selection pressures force the virus population to move in nearly a straight 
line in antigenic space ( Figure T^ ). Across replicate simulations, 94% of the variance of 



antigenic phenotype can be explained by a single dimension of variation. This mirrors the 
empirical results showing a largely linear antigenic map for H3N2 influenza isolates from 
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1968 to 2003 [9j. Because of the primarily one-dimensional movement, antigenic distance 
from the original phenotype increases nearly linearly with time (Figure "Sl^). Anti genic 
evolution occurs in a punctuated fashion; periods of relative stasis are interspersed with 
more rapid antigenic change ( Figure Sl^ ). Antigenic and epidemiological dynamics show 
a fundamental linkage, so that large jumps of antigenic phenotype result in increased rates 



of infection (Figure 1, Figure S2). In general, evolution via many smaller mutations results 



in a smoother antigenic map and less variation in yearly epidemics (Figure S3), while 



evolution via rare mutations of large effect exhibits a more clustered antigenic map and 



wider variation in seasonal incidence (Figure S4). 



The genealogical tree connecting the evolving virus population appears characteristically 
sparse with pronounced trunk and short-lived side branches ( Figure This tree shape 
is reflected in low levels of standing diversity; across replicates, an average of 5.68 years 
of evolution separate two randomly sampled viruses in the population. This level of di- 
versity matches what is observed in phylogenies of influenza A (H3N2) [T2|. A spindly 
genealogical tree is indicative of population turnover, wherein novel antigenic phenotypes 
continually replace more primitive 'spent' phenotypes, purging their genealogical diversity. 
In general, natural selection reduces effective population size and genealogical diversity [13]. 
Here, by comparing mutations occurring on trunk branches vs. mutations occurring on side 



branches, we find evidence for pervasive positive selection for antigenic change (Table 1) 



Trunk mutations tend to push antigenic phenotype forward along the line of primary anti- 



genic variation (Figure S5). We find a roughly linear relationship between the antigenic 



effect of a mutation and the likelihood that this mutation becomes incorporated into the 
trunk (Figure S6). Additionally, we find that trunk mutations occur at strikingly regular 



intervals, with less variation of waiting times than expected under a simple random process 



(Figure S7). There is a relative scarcity of mutation events occurring in intervals under 1 



year and a relative excess of a mutation events occurring in 2-3 year intervals (Figure S7). 



The genealogical tree also contains detailed information on the history of migration between 
regions. We find that, consistent with empirical estimates [HI US], the trunk resides pri- 



marily within the tropics, where seasonal dynamics are less prevalent (Figure 24). Across 



replicate simulations, we observe 72% of the trunk's history within the tropics and 28% 
within temperate regions. With symmetric host contact rates and equal host population 
sizes, and without seasonal forcing, we would expect trunk proportions of one third for each 
region. We calculated rates of migration based on observed event counts across replicate 
simulations, separating region-specific rates on side branches from region-specific rates on 
trunk branches. We find that migration patterns on side branches are close to symmet- 
ric, with similar rates between all regions, while migration patterns on trunk branches are 
highly asymmetric, with high rates of movement between temperate regions and from tem- 
perate regions into the tropics ( Figure 2^ ). Extrapolating from these rates, we arrive at 
an expected stationary distribution of 76% tropics and 24% temperate regions, in line with 
the observed residency patterns of the trunk. It may at first seem counter-intuitive to see 
higher rates of movement from the temperate regions into the tropics along trunk branches, 
but it makes sense when thought of in terms of conditional probability. Only those lineages 
that migrate into the tropics or those lineages which rapidly migrate between the north 
and south have a chance at becoming the trunk lineage, while lineages that remain within 
the temperate regions are doomed to extinction. 
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These findings suggest that persistence and migration are fundamentaUy connected and 
have important implications for future phylogeographic analyses. Russell et al. [Tl] em- 
phasize a source-sink model of movement of the HA protein of influenza A (H3N2) based 
on their finding of a trunk lineage residing within China and the Southeast Asian tropics. 
Whereas, Bedford et al. |15j emphasize a global metapopulation model based on phylo- 
genetic inference of migration rates across the entire tree. Our results suggest that both 
scenarios are simultaneously possible; side branches may be highly volatile moving rapidly 
and symmetrically between regions, while the trunk lineage may be more stable remaining 
within a region (or within a highly connected network of regions) that has more persistent 
transmission. In light of these results, we suggest that future work on the phylogeography 
of influenza take into account trunk vs. side branch differences in migration patterns. 



Discussion 



Correspondence between model and data 



Although multiple epidemiological/evolutionary mechanisms have been proposed to ex- 
plain the restricted genetic diversity and rapid population turnover of influenza A (H3N2) 
[H El m [7] , our results show that a simple model coupling antigenic and genealogical evo- 
lution exhibits broad explanatory power. We find a strong correspondence between the 



antigenic and genealogical patterns generated by our model ( Figure 1 ) and patterns of 



genetic and antigenic evolution exhibited by influenza A (H3N2) |3il2|- Our model sug- 
gests that punctuated antigenic evolution need only be explained by a lack of mutational 
opportunity and predicts that more detailed classification of influenza strains will support 
a relatively small number of predominant phenotypes ( Figure 1^ ). We suggest that a large 
proportion of intra-cluster variation in the observed antigenic map is due to experimental 
noise, rather than each strain possessing a unique antigenic location. Additionally, our 
model accurately predicts the contrasting dynamics of other types/subtypes of influenza. 
We find that lowering mutation size/effect or lowering intrinsic Rq results in decreased 
incidence, slower antigenic movement and greater genealogical diversity, all distinguishing 



characteristics of HlNl influenza and influenza B (Figure S8, Figure S9) 



In our model, when antigenic phenotype remains static, there may be multiple consecutive 



seasons without appreciable incidence ( Figure 1 A.) , a pattern apparently absent from H3N2 
influenza 1161 



We suggest that any model exhibiting punctuated evolution broadly consis- 
tent with the punctuated change seen in the experimental antigenic map will show similar 
patterns of incidence. We can 'fix' the incidence patterns, but at the cost of too smooth an 



antigenic map ( Figure S3 ) . Evolutionary patterns of the neuraminidase (NA) protein may 
provide an explanation. Epitopes in the HA and NA proteins are jointly responsible for 
determining antigenicity [2] , and it is now clear that levels of adaptive evolution are similar 
between HA and NA [T7]. Thus, changes in NA may be driving incidence patterns as well, 
resulting in an observed timeseries of incidence partially divorced from the antigenic map 
of HA. 



It remains a central question as to the extent that short-lived strain-transcending immu- 
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nity is responsible for influenza's limited diversity and spindly genealogical tree [U [5] . Our 
findings suggest a possible resolution. Although lacking short-lived immunity, our model 
shows a detailed correspondence to both the antigenic map and genealogical tree of H3N2 
influenza. If an antigenic map were to show a deep bifurcation, where two viral lineages 
move in different antigenic directions, then we would expect the same bifurcation to be 
evident in the genealogical tree. Short-lived strain-transcending immunity provides a mech- 
anism by which lineages may diverge in antigenic phenotype, but still show epidemiological 
interference. This mechanism would explain a situation where bifurcations emerge in the 
antigenic map, but competition results in the extinction of divergent antigenic lineages. 
The empirical antigenic map |9J suggests that this is not the case; one cluster leads to 
another cluster in orderly succession and there is never competition between antigenically 
distant clusters. This supports the hypothesis that antigenic evolution is primarily limited 
by a lack of mutational availability. This is not to say that short-lived strain-transcending 
immunity is not present; observed interference between subtypes [U [18] and evolution at 
CTL epitopes [19] provides substantial evidence for its existence. Instead, we suggest that 
short-lived strain-transcending immunity does not automatically generate antigenic maps 
and genealogical trees consistent with empirical evidence. 



Linear antigenic movement 

It would seem possible for one viral lineage to move in one antigenic direction, while an- 
other lineage moves tangentially, eventually resulting in two non-interacting viral lineages. 
Instead, we find that only movement in a single antigenic direction is favored. The origins 
of this pattern can be seen in the interaction between virus evolution and host immunity 



(Figure 3). As the virus population evolves forward it leaves a wake of immunity in the 
host population, and evolution away from this immunity results in the canalization of the 
antigenic phenotype; mutations that continue along the line of primary antigenic variation 
will show a transmission advantage compared to more tangential mutations. 

Following the work of Smith et al. [9], it remained an open question of why a two- 
dimensional map should explain the antigenic variation of H3N2 influenza. Although the 
authors astutely speculated that "there is a selective advantage for clusters that move away 
linearly from previous clusters as they most effectively escape existing population-level im- 
munity, and this is a plausible explanation for the somewhat linear antigenic evolution in 
regions of the antigenic map." This hypothesis remained to be tested. Here, we show from 
a simple model of epidemiology and evolution that a linear trajectory of antigenic evolution 
dynamically emerges due to basic selective pressures. This result simultaneously explains 
the linear pattern of antigenic drift [9| and the characteristically spindly genealogical tree 
^ exhibited by influenza A (H3N2). 

To consider to what extent these results were contingent on the dimensionality of the un- 
derlying antigenic model, we further implemented our model in a 10-dimensional antigenic 
space. Here, mutations occur as 10-spheres, but the distance moved by a mutation is the 
same as in the previous two-dimensional formulation. We arrive at nearly the same results 
with this model; principal components analysis shows that the first and second dimensions 



of variation account for 87% and 7%, respectively, of the total variance (Figure SIO). Thus 



6 



our model predicts that future work probing mutational effects will support an underly- 
ing high-dimensional antigenic space, even though a two-dimensional map is sufficient to 
explain observed antigenic relationships among strains. 



Winding back the tape 



It seems clear that, in our model, selection reduces the degrees of freedom of antigenic 
evolution. In light of this, we wanted to examine the degree of stochasticity in replicated 
evolutionary trajectories, and thereby test what happens when we "wind back the tape" 
|20| on the evolution of the virus. We ran 100 replicate simulations, each starting from 
the endpoint of the original 40-year simulation (Figure 4, Figure 5). Initially, we find 



a great detail of repeatability; during the first year of evolution, every replicate virus 
population undergoes a similar antigenic transition (Figure 4), resulting in a repeatable 



peak in northern hemisphere incidence (Figure 5). After three years, repeatability has 
mostly disappeared, with antigenic phenotype and incidence appearing highly variable 
across replicates (Figure 4 Figure 5). The 1-2 year timescale of repeatability can be 



explained by the presence of standing antigenic variation. In the initial virus population. 



there are several novel antigenic variants present at low frequency (Figure 3), one of which, 
without fail, comes to predominate the virus population. 

We see that the initial evolutionary trajectory, during which time standing variation plays 
out, is highly repeatable, and thus predictable given enough information and the right 
methods of analysis. However, prediction of longer-term evolutionary scenarios will neces- 
sarily be difficult or impossible except in a vague sense. Through careful surveillance efforts 
and genetic and antigenic characterization of influenza strains, the World Health Organi- 
zation makes twice-yearly vaccine strain recommendations [21]. It should be possible to 
combine these sorts of modeling approaches with surveillance data to gauge the likelihood 
that a sampled variant will spread through the population. 

Recent work on empirical fitness landscapes has shown that natural selection follows few 
mutational paths [22]. The spindly genealogical tree and the almost linear serial replace- 
ment of influenza strains has remained a puzzling phenomenon. We suggest that the 
evolutionary and epidemiological dynamics displayed by the influenza virus may simply 
be explained as an outgrowth of selection to avoid host immunity. Natural selection can 
only 'see' one step ahead, and so sacriflces long-term gains for short-term advantages. The 
result is a canalized evolutionary trajectory lacking antigenic diversiflcation. 



Materials and Metliods 

Transmission model 

To characterize the joint epidemiological, genealogical, antigenic and spatial patterns of 
influenza, we implemented a large-scale individual-based model. This model consists of 
daily time steps, in which the states of hosts and viruses are updated. Hosts may be born, 
may die, may contact other hosts allowing viral transmission, or may recover from infection. 
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Viruses may mutate in antigenic phenotype. Each simulation ran for 40 years of model 
time. 

Hosts in this model are divided between three regions: North, South and Tropics. There are 
30 million hosts within each of the three regions, giving = 9 x 10^ hosts. Host population 
size remains fixed at this number, but vital dynamics cause births and deaths of hosts at a 
rate of 1/30 years = 9.1 x 10~^ per host per day. Within each region, transmission proceeds 
through mass- action with contacts between hosts occurring at a rate of (3 = 0.36 per host 
per day. Regional transmission rates in temperate regions vary according to sinusoidal 
seasonal forcing with amplitude e = 0.15 and opposite phase in the North and in the 
South. Transmission rate does not vary over time in the Tropics. Transmission between 
region i and region j occurs at rate m/3i, where m = 0.001 and is the same between each 
pair of regions and /3j is the within-region contact rate. Hosts recover from infection at 
rate v = 0.2 per host per day, so that Rq in a naive host population is 1.8. There is no 
super-infection in the model. 

Each virus possesses an antigenic phenotype, represented as a location in Euclidean space. 
Here, we primarily use a two-dimensional antigenic location. After recovery, a host 're- 
members' the phenotype of its infecting virus as part of its immune history. When a 
contact event occurs and a virus attempts to infect a host, the Euclidean distance from 
infecting phenotype (p^ is calculated to each of the phenotypes in the host immune history 
, . . . , (j)h„ ■ Here, one unit of antigenic distance is designed to correspond to a twofold 
dilution of antiserum in a hemagglutination inhibition (HI) assay [9] . The probability that 
infection occurs after exposure is proportional to the distance d to the closest phenotype in 
the host immune history. Risk of infection follows the form p = minjd s, 1}, where s = 0.07. 
Cross-immunity a equals 1 — p. The initial host population begins with enough immunity 
to slow down the initial virus upswing and place the dynamics closer to their equilibrium 
state; initial R was 1.28. 

Our model follows Gog and Grenfell |23j in representing antigenic distance as distance 
between points in a geometric space. By forcing one-dimension to directly modulate /3, 
Gog and Grenfell find an orderly replacement of strains. Kryazhimskiy et al. 124] use 
a two-dimensional strain-space, but enforce a cross-immunity kernel that directly favors 
moving along a diagonal line away from previous strains. Our model does not 'build in' 
the one-dimensional direction of antigenic drift, which instead emerges dynamically from 
competition among strains. 

The initial virus population consisted of 10 infections each with the identical antigenic 
phenotype of {0,0}. Over time viruses evolve in antigenic phenotype. Each day there is 
a chance /x = 10~^ that an infection mutates to a new phenotype. This mutation rate 
represents a phenotypic rate, rather than genetic mutation rate, and can be thought of 
as arising from multiple genetic sources. When a mutation occurs, the virus's phenotype 
is moved in a completely random direction ~ Uniform(0, 360) degrees. Mutation size 
is sampled from the distribution ~ Gamma(a,/3), where a and /? are chosen to give a 
mean mutation size of 0.6 units and a standard deviation of 0.4 units. This distribution 
is parameterized so that mutation usually has little effect on antigenic phenotype, but 
occasionally has a large effect. This is similar to the neutral networks implemented by 
Koelle et al. [6j , wherein most amino acid changes result in little decrease to cross-immunity 
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between strains, but some changes result in large jumps in cross-immunity. 



Model output 

Daily incidence and prevalence are recorded for each region. During the course of the simu- 
lation, samples of current infections are taken from the evolving virus population at a rate 
proportional to prevalence. Each viral infection is assigned a unique ID, and in addition, 
infections have their phenotypes, locations and dates of infection recorded. In this model, 
viruses lack sequences and so standard phylogenetic inference of the evolutionary relation- 
ships among strains is impossible. Instead, the viral genealogy is directly recorded. This is 
made possible by tracking transmission events connecting infections during the simulation; 
infections record the ID of their 'parent' infection. Proceeding from a sample of infections, 
their genealogical history can be reconstructed by following consecutive links to parental 
infections. During this procedure, lineages coalesce to the ancestral lineages shared by 
the sampled infections, eventually arriving at the initial infection introduced at the be- 
ginning of the simulation. Commonly, phylodynamic simulations generate sequences that 
are subsequently analyzed with phylogenetic software to produce an estimated genealogy 
[UEIISS]. This step of phylogenetic inference is imperfect and computationally intensive, 
and by side-stepping phylogenetic reconstruction we arrive at genealogies quickly and ac- 
curately. Other authors have implemented similar tracking of infection trees |26| [27] . This 
genealogy-centric approach makes many otherwise difficult calculations transparent, such 



as calculating lineage-specific region-specific migration rates ( Figure 2 ) and lineage-specific 



mutation effects (Figure S5) 



Infections are sampled at a rate designed to give approximately 6000 samples over the 
course of the simulation, with genealogies constructed from a subsample of approximately 



300 samples. The results presented in Figure 1 represent a single representative model 



output; one hundred replicate simulations were conducted to arrive at statistical estimates. 



Parameter selection and sensitivity analysis 

Estimating what the basic reproductive number Rq for seasonal infiuenza would be in a 
naive population is notoriously difficult. Season-to-season estimates of effective reproduc- 
tive number R for the USA and France gathered from mortality timeseries display an 
interquartile range of 0.9-1.8 [2H]. Geographic spread within the USA suggests an average 
seasonal R of 1.35 |29j . These estimates of R will be lower than the Rq of influenza due to 
the effects of human immunity. We assumed Rq of 1.8, consistent with the upper range of 
seasonal estimates. Duration of infection was chosen based on patterns of viral shedding 
shown during challenge studies |30| . The linear form of the risk of infection and its increase 
as a function of antigenic distance s was chosen as 0.07 based on experimental work on 
equine influenza [31j and from studies of vaccine effectiveness (32j. Between-region contact 
rate m was chosen to yield a trunk lineage that resides predominantly in the tropics. With 
much higher rates of mixing, the trunk lineage ceases to show a preference the tropics, and 
with much lower rates of mixing, particular seasons in the north and the south will often 
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be skipped. The amplitude of seasonal forcing e was chosen to be just large enough to get 
consistent fade-outs in the summer months and is consistent with empirical estimates [33]. 

Mutational parameters were based, in part, on model behavior. We assumed 10 amino 
acid sites involved in antigenicity, each mutating at a rate of 10~^ [12] to give a phenotypic 
mutation rate /i = 10~^ per infection per day. We chose mutational effect parameters 
(mean = 0.6, sd = 0.4) that would give suitably fast rates of antigenic evolution corre- 
sponding to approximately 1.2 units of antigenic change per year, while simultaneously 
giving clustered patterns of antigenic evolution [9J. Similar outcomes are possible under 
a variety of parameterizations. If mutations are more common = 3 x 10^^) and show 
less variation in effect size (mean = 0.6, sd = 0.2), then antigenic drift occurs in a more 
continuous fashion, resulting in less variation in seasonal incidence and a smoother distri- 
bution of antigenic phenotypes (Figure S3). If mutations are less common (/x = 5 x 10~^) 



and show more variance in effect (mean = 0.7, sd = 0.5), then antigenic change occurs in 



a more punctuated fashion (Figure S4). Basic reproductive number Rq can be traded off 
with mutational parameters to some extent. Less mutational input and higher Rq will give 
similar patterns of antigenic drift and seasonal incidence. Similarly, Kucharski and Gog 
|34j find that increasing Rq results in increased rates of emergence of antigenically novel 
strains. 

In 20 out of the 100 replicate simulations, we observed a major bifurcation of antigenic 
phenotype and a consequent increase in incidence and genealogical diversity. These sim- 
ulations were removed from the analysis. Similar to Koelle et al. [35], we assume that 
although the historical evolution of H3N2 influenza followed the path of a single lineage, 
it could have included a major bifurcation. Further work in these directions will help to 
determine the likelihoods of single lineage vs. bifurcating scenarios. 



Antigenic map 

Antigenic phenotypes are modeled as discrete entities on the Euclidean plane; multiple 
samples have the same antigenic location. However, in the empirical antigenic map of 
influenza A (II3N2), each strain appears in a unique location [9]. We would argue that 
some of this pattern comes from experimental noise. Indeed, Smith et al. |9] find that 
observed measurements and measurements predicted from the map differ by an average 
of 0.83 antigenic units with a standard deviation of 0.67 antigenic units. We take this 
as a proxy for experimental noise and add jitter to each sampled antigenic phenotype by 
moving it in a random direction for an exponentially distributed distance with mean of 0.53 
antigenic units. If two samples with the same underlying antigenic phenotype are jittered 
in this fashion, the distance between them averages 0.83 antigenic units with a standard 
deviation of 0.64 units. 

We added noise to each of the 5943 sampled viruses in this fashion resulting in an approx- 
imated antigenic map ( Figure Ip ). Virus samples were then clustering following standard 



clustering algorithms. We tried clustering by the A;- means algorithm and also agglomer- 
ative hierarchical clustering with a variety of linkage criterion. We found that clustering 
by Ward's criterion consistently outperformed other methods, when measured in terms of 
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within-cluster and bctwccn-cluster variances. However, the exact clustering algorithm had 
little effect on our overall results. 
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Tables 



Table 1. Rates of mutation and phenotypic change on trunk and side branches and 
mutational expectation. 





Baseline 


Side branch 


Trunk 


Trunk / side branch 


Mutation size (AG units) 


0.60 


0.79 


1.58 


1.99X 


Mutation rate (mut per year) 


0.04 


0.06 


0.81 


13.23X 


Antigenic flux (AG units per year) 


0.02 


0.05 


1.27 


26.25X 
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Figures 




Figure 1. Simulation results showing epidemiological, antigenic and genealogical dy- 
namics. (A) Weekly timeseries of incidence of viral infection in north and tropics regions. (B) 
Two-dimensional antigenic phenotypes of 5943 viruses sampled over the course of the simulation. 
Each discrete virus phenotype is shown as a bubble, with bubble area proportional to the number 
of times this phenotype was sampled. (C) Genealogical tree depicting the infection history of 376 
samples from the virus population. Parent/offspring relationships were tracked over the course of 
the simulation, giving a direct observation of the genealogy rather than a phylogenetic inference. 
(D) Antigenic map depicting phenotypes of 5943 viruses sampled over the course of the simulation. 
To construct the map, noise was added to each sample and the resulting observations grouped into 
11 clusters and colored accordingly. Grid lines show single units of antigenic distance. Cluster 
assignments were used to color all panels in a consistent fashion. 



15 




Figure 2. Patterns of geographic movement of virus lineages. (A) Evolutionary relation- 
ships among 376 viruses sampled evenly through time colored by geographic location. Lineages 
residing in the north (N), south (S) and tropics (T) are colored yellow, red and blue respectively. 
(B) Observed migration rates between regions on side branch lineages (left) and on trunk lineages 
(right). Arrows denote movement of lineages and arrow width is proportional to migration rate. 
Circle area is proportional to the expected stationary frequency of a region given the observed 
migration rates. In both cases, migration rates are calculated across 80 replicate simulations. 



16 




26 28 30 32 34 36 38 40 
Antigenic dimen 1 



Figure 3. Host immunity and antigenic history of the virus population. Contour lines 
represent the state of host immunity at the end of the 40-year simulation. They show the mean 
risk of infection (as a percentage) after a random host in the population encounters a virus bearing 
a particular antigenic phenotype. Contour lines are spaced in intervals of 2.5%. Bubbles represent 
a sample of antigenic phenotypes present at the end of the 40-year simulation. The area of each 
bubble is proportional to the number of samples with this phenotype. Lines leading into these 
bubbles show past antigenic history. The current phenotypes rapidly coalesce to a trunk phenotype. 
The movement of the virus population from the left to the center of the figure can be seen from the 
antigenic history. At the end of the simulation several virus phenotypes exist with similar antigenic 
locations; all of these phenotypes lie significantly ahead of the peak of host immunity. 
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Figure 4. Antigenic phenotypes over the course of 4 years of evolution across 100 
replicate simulations starting from identical initial conditions. Replicate simulations were 
initialized with the end state of the 40-year simulation shown in Figure 1 Each panel shows an 



additional year of evolution, with black points representing the mean antigenic phenotypes of the 
100 replicate simulations and gray lines representing the history of each mean antigenic phenotype. 
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Figure 5. Timeseries of incidence across 100 replicate simulations with identical initial 
conditions. Panels show incidence in the North, Tropics and South regions over the course of 6 
years. Solid black lines represent the median weekly incidence across the 100 replicate simulations, 
while gray intervals represent the interquartile range across simulations. There is little variability 
for the first year of replicate simulations. Replicate simulations were initialized with the end state 



of the 40-year simulation shown in Figure 1 
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Figure SI. Antigenic evolution over the course of the 40-year simulation. (A) Proportion 
of virus population comprised of each antigenic cluster through time. (B) Antigenic distance from 
initial phenotype {x = 0, y = 0) for each of 5943 virus samples relative to time of virus sampling. 
Viruses were sampled at a constant rate proportional to prevalence and coloring was determined 



from the antigenic map in Figure 1 3 
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Antigenic drift between years f and f+1 

Figure S2. Correlation between antigenic drift and attack rate. Antigenic drift is measured 
as the distance between tlie centroid of plienotypes of one year and the centroid of phenotypes of 
the following year. Measurements were taken across 80 replicate simulations. Individual pairs of 
measurements are shown as gray points and a locally-linear regression (LOESS) is shown as a black 
dashed line. 
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Figure S3. Simulation results showing epidemiological, antigenic and genealogical dy- 
namics for 'smoother' mutation model. (A) Weekly timeseries of incidence of viral infection 
in north and tropics regions. (B) Antigenic map depicting phenotypes of viruses sampled over the 
course of the simulation. Grid lines show single units of antigenic distance. (C) Genealogical tree 
depicting the infection history of samples from the virus population. Cluster assignments were used 
to color panels (A), (B) and (C) in a consistent fashion. Alternative mutational parameters are 
/i = 3 X 10~^, mean mutation size of 0.6 units and standard deviation of mutation size of 0.2 units. 
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Figure S4. Simulation results showing epidemiological, antigenic and genealogical dy- 
namics for 'rougher' mutation model. (A) Weekly timeseries of incidence of viral infection 
in north and tropics regions. (B) Antigenic map depicting phenotypes of viruses sampled over the 
course of the simulation. Grid lines show single units of antigenic distance. (C) Genealogical tree 
depicting the infection history of samples from the virus population. Cluster assignments were used 
to color panels (A), (B) and (C) in a consistent fashion. Alternative mutational parameters are 
/i = 5 X 10~^, mean mutation size of 0.7 units and standard deviation of mutation size of 0.5 units. 
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Figure S5. Mutation spectrum in two-dimensional antigenic space of side branch mu- 
tations and trunk mutations. (A) Histogram of mutation effects along the axis of primary 

antigenic variation across 80 replicate simulations. The left panel shows the distribution of effects 
of side branch mutations and the right panel shows the distribution of effects of trunk mutations. 
(B) Smoothed two-dimensional histogram of mutation effects along the primary and secondary axes 
of antigenic variation across 80 replicate simulations. Histograms were constructed from 21,405 side 
branch mutations and 1584 trunk mutations. 
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Figure S6. Relationship between a mutation's phenotypic effect and its likelihood of 
being part of the trunk. The a;-axis represents the effect of a mutation along the line of primary 

antigenic variation, and the j/-axis represents the probability that the mutation is part of the trunk. 
Mutations of large effect are increasingly rare, but when they do occur are increasingly likely to be 
part of the trunk. 
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Figure S7. Observed vs. expected distributions of waiting times between phenotypic 
mutations along genealogy trunk. (A) Histogram bins show the observed distribution of wait- 
ing times in years across 80 replicate simulations representing 1584 mutations. The mean of this 
distribution is 1.76 years. The dashed line shows the Poisson process expectation of exponentially 
distributed waiting times. (B) The density distribution of waiting times is transformed into a haz- 
ard function, representing the rate of trunk mutation after a specific waiting time. The dashed line 
shows the memoryless hazard function of the Poisson process expectation. 
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Figure S8. Simulation results showing epidemiological, antigenic and genealogical dy- 
namics with weaker mutation. (A) Weekly timeseries of incidence of viral infection in north 
and tropics regions. (B) Antigenic map depicting phenotypes of viruses sampled over the course of 
the simulation. Grid lines show single units of antigenic distance. (C) Genealogical tree depicting 
the infection history of samples from the virus population. Cluster assignments were used to color 
panels (A), (B) and (C) in a consistent fashion. Here, /i = 5 x 10~^, mean mutation size is 0.42 
units and standard deviation of mutation size is 0.28 units. 
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Figure S9. Simulation results showing epidemiological, antigenic and genealogical dy- 
namics with lower intrinsic i?o- (A) Weekly timeseries of incidence of viral infection in north 
and tropics regions. (B) Antigenic map depicting phenotypes of viruses sampled over the course of 
the simulation. Grid lines show single units of antigenic distance. (C) Genealogical tree depicting 
the infection history of samples from the virus population. Cluster assignments were used to color 
panels (A), (B) and (C) in a consistent fashion. Here, (3 = 0.3, giving i?o = 1-5. 
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Figure SIO. Principal components of antigenic variation under a 10-sphere mutation 
model. Each panel shows 5991 samples of antigenic phenotype over the course of a 40-year simu- 
lation. Each phenotype is represented as a bubble, with bubble area proportional to the number of 
samples with this phenotype. Bubbles are colored based on clustering the 10-dimensional antigenic 
phenotypes. The original 10-dimensional space was rotated using principal components analysis to 
give orthogonal axes in the order of their contribution to antigenic variation. Each panel shows 
a two-dimensional slice of the this rotated space. Principal components 7-10 were left out of the 
figure for clarity. 
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